Strongly-coupled Josephson junction array for simulation of frustrated 

one-dimensional spin models 
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We study the capacitance-coupled Josephson junction array beyond the small-coupling limit. We 
find that, when the scale of the system is large, its Hamiltonian can be obtained without the small- 
coupling approximation and the system can be used to simulate strongly frustrated one-dimensional 
Ising spin problems. To engineer the system Hamiltonian for an ideal theoretical model, we apply 
a dynamical decoupling technique to eliminate undesirable couplings in the system. Using a 6-site 
junction array as an example, we numerically evaluate the system to show that it exhibits important 
characteristics of the frustrated spin model. 

PACS numbers: 03.67.Ac, 75.10.Jm, 85.25.Cp 



I. INTRODUCTION 

As an important application of quantum information 
science, quantum simulation of difficult physics problems 
has received much attention in recent years. Theoreti- 
cally, there have been many proposals of quantum sim- 
ulators based on various physical systems [iHal- Experi- 
mentally, simulations of some important physics models 
have been demonstrated (6l-[Tl|. 

Quantum simulation is most valuable for studying 
strongly correlated problems since there are no generally- 
applicable theoretical methods to solve them. The strong 
interactions involved in these problems usually translate 
into strong couplings between entities in a quantum sim- 
ulator. This requirement of strong couplings often poses 
a challenge for the design and construction of a quan- 
tum simulator, since it can be difficult to engineer such 
couplings in a simulation system. Even when strong cou- 
plings are available, it can still be nontrivial to tailor the 
couplings in an well-controlled manner that is required 
for the problems to be simulated. As one such example, 
it is usually difficult to obtain the system Hamiltonian for 
a Josephson junction array when the couplings between 
the junctions are strong. Consequently, Josephson de- 
vice based quantum simulation systems often operate in 
the small coupling limit in which the system Hamiltonian 
can be obtained by treating the coupling as perturbation 

In order to go beyond the small-coupling limit and con- 
struct a system useful for simulating strongly correlated 
physics, in this paper we investigate a one-dimensional 
Josephson junction array which is coupled by large ca- 
pacitances that cannot be treated perturbatively. Inter- 
estingly, the system Hamiltonian can be obtained exactly 
without the small coupling approximation. It is found 
that, in the large coupling limit, the interaction strength 



between next nearest neighbors can become comparable 
with that between the nearest neighbors. Because of this, 
we can use the system to study the important problem 
of one-dimensional frustrated spin models whose phase 
diagrams and properties have not been completely re- 
solved (TB - fTol . In order to control the system Hamil- 
tonian to match that of the ideal theoretical model, we 
use a dynamical decoupling technique to suppress inter- 
actions between neighbors that are three site locations or 
farther apart. 



II. SYSTEM HAMILTONIAN OF THE 
JOSEPHSON-JUNCTION ARRAY 

The system we study is an A'^-site Josephson junction 
array as shown in Fig. [T] Each site consists of a charge 
island biased by a voltage source V^. through a gate ca- 
pacitance Cg., where i = 1,...,A^. The charge on the 
island can tunnel through a SQUID device whose total 
capacitance is Cj and whose effective Josephson energy 
Ej can be adjusted by a flux bias. Adjacent charge is- 
lands, as well as those at the ends of the array, are cou- 
pled by a capacitance Cc- We take the average phase 
ifi of the SQUID on site i as the generalized coordinates 
for the system. Its rate of change is determined by the 
voltage Vj across the Josephson junction according to the 
Josephson relation Vj = {h/2e)(pi [20]. The charging en- 
ergy T of the capacitances and the Josephson energy V 
of the junctions are 
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From the Lagrangian £ = T — V of the system, we can 
derive the generalized momentum which is related to the 
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charge number on the islands 

Pi = dC/dipi = -hui 



(3) 



Denoting the charge number on the ith island rii and 
the bias charge number Ugi = Cg-Vg-/2e, we can write 
the system's equation of motion 



(4) 



where tp = 
'lg2, ■■■,ni 



[(pi,(p2,--,(Pi,-;iPN\ 

Igi, TIN 



[ni - ngi,n2 



n„2, rii — n„i, njv — ng^]^ , and M is the matrix 
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2Cc is the total capacitance of 



where Cs = Cj 
each island. 

The Hamiltonian of this Josephson junction array is: 

i 

1 , h 



-{2e)^n^M-^n + V. 



(6) 



In order to obtain the system Hamiltonian, the inverse 
matrix of M must be calculated. Since it is nontrivial 
to exactly solve for M~^, most previous studies [l2l. [l3l| 
have assumed a small coupling capacitance Cc so the sys- 
tem Hamiltonian can be obtained by treating the cou- 
pling as perturbation. Unfortunately, this precludes the 
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FIG. 1: (Color online) The capacitance-coupled Josephson 
junction array. Cc is the coupling capacitance between neigh- 
boring junctions. Vgi is the bias voltage and Cgi is the bias 
capacitance of the Josephson junction at site i. The total 
capacitance and Josephson energy are Cji and £j. 



system from being used to simulate strongly-correlated 
problems in which strong couplings are required. Here, 
we assume that Cc is not necessarily small and try to 
solve for exactly. 

Considering the translational symmetry of the prob- 
lem, we see the inverse matrix must be in the form 
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, (7) 



Using some mathematic techniques for solving polyno- 
mials, we can calculate the values of the matrix elements 
in exactly: 



o,^-^{X^-'Ao + ^Bo), 



where 



Ao = 



A = 



l-2/3A + (l-f)^A^-i' 
(2^- A)A^-i 



1 - 2/3A 
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2Cc + Cgi + Cj 2 
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(11) 
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When N ^ 00, the above results simplify to ^ 
1/(1 - 2^A), Bo 0, and X'-^Aq/Cs- 

As can be seen in Eq. A characterizes the ratio 
between adjacent and non-adjacent interaction strengths 
in our system. According to Eqs. (|11|) and (fT2|). A 
is determined by the coupling capacitance Cc- In the 
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FIG. 2: (Color online) The relationship between A and /3. The 
dashed red line corresponds to A = 1/2. 



weak coupling limit Cc <C Cs, /3 <C 1 and A is nearly 
equal to j3. However, when the coupling is strong, A in- 
creases quickly, as shown in Fig. [2] In particular, when 
the coupling capacitance Cc dominates, A can approach 
1, and the Hamiltonian in Eq. (|6]) describes a deeply 
frustrated system with appreciable non-adjacent interac- 
tions. Meanwhile, according to the results of Ref. [2l| . 
there exist many close energy levels as /? approaches 1/2 
too closely, which can fail the two-level approximation. 
Therefore we only pay our attention to a large /3, but not 
closely approaching 1/2 in the rest of our paper. 

Under proper conditions, if we bias the charge islands 
at the vicinity of Ugi = 1/2, we can use the two-level 
approximation for the charge qubits with tii = Q and 
ni = 1 as the basis states. We can then write the system 
Hamiltonian in the following Pauli matrix representation 



where the spin up and down states represent the = 
and rii = 1 states, and the transverse magnetic field 
B~—8jC^/{2\e'^A[)). B can be adjusted by tuning the 
magnetic flux of the SQUIDs. Here we applied a canon- 
ical cr^ transformation on even sites for the convenience 
of our following discussion about phase diagram without 
changing the physics of the system. Usually, two-level ap- 
proximation works very well for a single Josephson charge 
island [23 - [23 |. In our system of Josephson junction array, 
more careful analysis is necessary. Appendix A provides 
a detailed discussion about the applicability of the two- 
level approximation in our system. 



III. SIMULATION OF THE ANNNI MODEL 
USING DYNAMICAL DECOUPLING 

Our circuit is useful for quantum simulation of frus- 
trated spin problems since it exhibits strong non-adjacent 
spin interactions in the large coupling limit. However, 
the Hamiltonian in Eq. (fT3|) does not correspond to an 
ideal theoretical model with limited-range interactions 
yet. We intend to further engineer it for quantum simu- 
lation of well known frustrated models. As an example, 
we show how to simulate the one dimensional axial next- 
nearest-neighbor Ising (ANNNI) model in external fields. 
Its Hamiltonian is given by 

i i i 

The ANNNI problem is an important model for studying 
frustrated physics due to competition between adjacent 
and non-adjacent neighbor interactions [25, ,26|. Despite 
years of research, its phase diagrams and physic al prop- 
erties have not been completely resolved |15l-[l9j. 

Comparing our circuit Hamiltonian in Eq. (I13|) and 
target Hamiltonian in Eq. (jl4p . we find that there are 
extra terms that describe interactions between spins that 
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FIG. 3: (Color online) The dynamical decoupling control 
scheme to eliminate interactions between spins that are 3 sites 
apart. The ellipses in groups of 3 along the vertical direction 
represent simultaneous n or 3n pulses applied to the qubits. 



are 3 or more sites apart. We eliminate these extra terms 
by using techniques of dynamical decoupling [13, ■ In 
this practice, we apply carefully designed sequences of 
fast short pulses to individual qubits to engineer a Hamil- 
tonian that can be very different than the original Hamil- 
tonian. 

In the following, we will demonstrate how to eliminate 
interactions between spins that are 3 sites apart for a 
spin chain with the Hamiltonian in Eq. (fT3)) . For clari- 
fication and ease of illustration, we discuss the details of 
our scheme on a 6-site spin chain with periodic bound- 
ary condition. The same technique applies in a spin chain 
with arbitrary length. 

When the number of qubits in the system is 6, the 
Hamiltonian in Eq. ([T3|) reads 

Hs = Ji{cFl(7l+a2(jl + dial + al(7l+alal + dial) 



(15) 



where Ji and J2 are interaction strengths between the 
nearest and next-nearest neighbors and J3 is that be- 
tween spins that are 3 sites apart. 

Our scheme to engineer the desired Hamiltonian 
involves applying rapid pulsed operations Rx{0) — 
exp{— (i0(T2;)/2} on individual qubits. When 6 — ir (37r) 
the corresponding unitary operation is Rxi"^) = — ?cr^ 
(i?2;(37r) = R\;{'k) — ia^). We have the following com- 
mutation relations: 



i?i(^)(T^i?,(^) = a% 



(16) 
(17) 

Rl{^)e-'"'Rx{^) = e-'[Rl(-)HRA-)]t_ (13) 

Our procedure contains four phases as shown in Fig. 
[3l Each phase is completed in an interval of 25t where 5t 
is a short time. At the beginning of each phase, the oper- 
ations i?F''-)(^) = i?^(^)i?^,(^)i?^(7r) (z, j, fc = 1, 2, 6) 
are applied to the qubits on sites and k simultane- 
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ously (see Fig. [3]). After a time interval of St, the con- 
jugate operations i?i^^'^'^'*(7r) are applied in the second 
half of the phase. The evolution of the system at the end 
of phase 1 is given by 

U, = e-^^»**i?t(l-2^3) ^(1.2,3) (^)_ (19) 

For a short time duration 6t, we have g^-'^HiSt^-iH^st _ 
g-j(Hi+//2)5t_^(9(^^2)^ To first order in 5t, the unitary op- 
erator Ui = exp{—iH'^-^-^2St), where the effective Hamil- 
tonian in phase 1 is 

= -^llo-Jo-f + cr|crf -t- cr|(T| -I- o-fcr|) 

6 

+J,{alal + al<J^,)~BY,af. (20) 

i=l 

From Eq. (|^D|) . we can see that all couplings between 
qubits that are 3 sites apart have been eliminated. How- 
ever, some terms that we want to keep, such as the 
nearest- neighbor coupling cr|(T| and next-nearest cou- 
pling cr|cr|, are also eliminated. In order to make up for 
this problem, in phase 2 and 3 we use the same technique 
but shift the target qubits one site a time as shown in Fig. 
[31 This gives us the following effective Hamiltonian for 
these two phases 

6 

+M<j^,al + al<jt)~BY,af; (21) 
1=1 

3 = -'l(CTif72 + <^3f^4 + C^4<^5 + ^e^^l) 

6 

^Ualal^alaD^BY^ot. (22) 
i=\ 

Obviously, the missing terms for nearest and next-nearest 
neighbor couplings in ^ in Eq. (|20p are compensated 
by remaining terms in Eqs. (|2ip and ^F2\ . Similarly, 
missing terms in H^^^ and H*^^^ are compensated. Nev- 
ertheless, some next-nearest neighbor coupling terms are 
still missing from the sum of H^-^-^ , H2^^ and H^^-^ , since 
in each phase 2/3 of the nearest- neighbor couplings re- 
main but only 1/3 of the next-nearest-neighbor couplings 
survive. In order to obtain all the nearest-neighbor and 
the next-nearest-neighbor couplings, we need a phase 4 
as shown in Fig. [3l By keeping the next-nearest-neighbor 
interactions and eliminating the nearest-neighbor inter- 
actions, it gives us the effective Hamiltonian 

= ^2(^10-3 +(T2Cr4 +C^3^5 +^4'^6 +f^5^1 +C^6f^2)) 
6 

-B^af. (23) 
The combined evolution of the system for the 4 phases 



is 

U = t/4f/3[/2f/i«e-^^4^"25tg-.H3=^^25tg-,H|^^25tg~.ff,^"25t 
A ^-iH^sfSSt (24) 

where the effective average Hamiltonian is 

HeSJ = -[Jiiulal + (J2<^1 + alal + alal + dial + alal) 

I 7 ( „z z I „z z I „z z I „z z I „z z I „z zw 
-^'h\P\<yz + Cr2^4+f^3'^5+<^4f^6+'^5'^l+'^6'^2)J 

-5E<- (25) 

i=\ 

This is exactly the ANNNI model with periodic boundary 
condition. 

Though we have used a 6-site chain to demonstrate 
how to eliminate couplings between qubits that are 3 
sites apart, it is obvious that, by applying the operations 
Rx''''''\'^) to groups of 3 qubits in the chain and shifting 
the target qubits by 1 site at a time in each phase, we can 
use the same technique to eliminate couplings between 
qubits 3 sites apart in an infinite-length spin chain. No- 
tice that in Eq. (|13p there are couplings between qubits 4 
or more sites apart. These couplings are weaker since the 
interaction strengths A-'^^ in Eq. p3)) decreases with the 
distance between qubits, but the error caused by them 
may still be unacceptable depending on the required ac- 
curacy of the simulation. By using a nested dynamical 
decoupling scheme, it can be shown that all couplings 
between qubits that are separated by 3 or more sites can 
be eliminated [1^. Therefore, given a required accuracy, 
we can in principle achieve the ANNNI Hamiltonian in 
Eq. (0. 

IV. PHASE DIAGRAM OF THE ANNNI 
MODEL 

Now that we can simulate the ANNNI model, we per- 
form some analysis on its phase diagram. When B = 
and A is small, the nearest-neighbor interactions dom- 
inate and the ground state is the ferromagnetic state 
• • • • • • )z (or I tt • • • tt • • • )z)- As A increases, 

the next-nearest-neighbor interactions become impor- 
tant. When A reaches some critical value, they become 
the dominating factor and the antiphase | ttiitt ' ' •)z 
(or I J^J-ttii ■ ■ ■ ) z ) which minimizes the next-nearest 
neighbor interaction energy becomes the ground state. 
In the limit of large transverse field B, the ground state 
will be the paramagnetic phase | tt ■ ■ ■ tt ■ ■ ' )a: to mini- 
mize the Zeeman energies. 

From the above analysis, we see that there should be 
a ferromagnetic phase, a paramagnetic phase and an an- 
tiphase in the ANNNI model. However, there could be 
more subtle regimes in the phase diagram. Studies have 
shown (inconclusively) that there could be a unique fioat- 
ing phase in the deeply frustrated regime [l^, [l^ . This 
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phase is characterized by the fact that the nth-neighbor 
spin-spin correlation function in the longitudinal direc- 
tion decays algebraically. The exact origin and range of 
the floating phase is still an open question and therefore 
a good subject for quantum simulation. Since the float- 
ing phase is located in the deeply frustrated regime, our 
circuit with strong couplings offers a good system for its 
simulation. 

Many numerical recipes such as finite-size scaling 
method [13 and the interface approach [l^ have been 
used to calculate the phase diagram of frustrated Ising 
model. We use the new method of time-evolving block 
decimation (TEBD) algorithm [s^, HH to calculate the 
ground state energy of ANNNI model. TEBD is an pow- 
erful algorithm to simulate quantum evolution process 
based on matrix product state and Trotter expansion 
[30I [3l| . By making the evolution time imaginary, we 
get the so called i-TEBD which can be used to deter- 
mine the ground state of a system efficiently. The TEBD 
method is based on the following matrix representation 
of a quantum state 

d d 

I*) = H • • • c^i-" 1*1) 1^") (26) 

where d is the number of local energy levels on every site. 
The coefficients 

X 

f.. . = \^ rl^l'i/^WrPI'^ /^[2]p[3]i3 ...-rNi„ (oT) 

ai,--- ,Qt^_i 

are defined with the help of n tensors {F^^l, • • • , r'"!} and 
n — 1 vectors ■ ■ ■ ,^'"~^'}, where x is the maximal 
number of two-party Schmidt decomposition coefficients. 
In practice, % does not need to be very large, because the 
Schmidt coefficients roughly decay exponentially with a. 
Any single-site operation or adjacent-site joined opera- 
tion on the state can be achieved by updating the corre- 
sponding tensors and vectors. Here, we use the second 
order Trotter expansion for i-TEBD (see reference [3l|). 
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FIG. 5: (Color online) (a) Fine structures and small dips in 
the dEg/dB^ curve with the parameters A = 0.7 and = 60. 
(b) The phase diagram of the ANNNI model, (c) The spin- 
spin correlation function (d) versus the spin separation d at 
point A (A = 0.7 and B — 0.2) in the phase diagram, (d) 
c^ld) at point B (A = 0.7 and B = 0.3) in the phase diagram, 
(e) Cs(ci) at point C (A = 0.7 and B = 0.6) in the phase 
diagram. 
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FIG. 4: (Color online) The finite-size scaling of second order 
derivative of ground state energy with A = (solid lines) and 
A = 0.2 (dashed lines) 

Using the TEBD algorithm, we calculate the ground 
state energy of ANNNI model in the external field B. 



How the ground state energy changes with the parame- 
ters in the Hamiltonian is of great significance because it 
gives us important clues on the quantum phase transition 
points. Considering this, we plot the second derivative 
of the ground state energy d^E/dB^ in Fig. |4l for differ- 
ent coupling strength A and system size N . Notice that 
there are dips in the curves in Fig. |4l Their positions 
nearly do not change with the system size N when N is 
large enough. These dips are where quantum phase tran- 
sitions occur, and we can read from their positions the 
corresponding critical field strengths at the phase transi- 
tion points. These critical parameter values allow us to 
construct the system's phase diagram which is shown in 
Fig. 5(b). The phase diagram is consistent with earlier 
results [1^ obtained using different numerical methods. 

Interestingly, as shown in Fig. 5(a) we find that in the 
strongly frustrated regime there can be a segment in the 
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curve of (PEg/dB"^ where there are muhiple extra sniall 
dips in addition to the main dip. This region is labeled 
by the green area in the phase diagram in Fig. 5(b). It 
is roughly at the location of the floating phase obtained 
in earlier work 19]. To further clarify the characteristics 
of the system in this region, we calculate the spin-spin 
correlation function 



cs{d)-- 



'N/2+l^N/2+l+d/ YN/2+1 /YN/2+l+d 



(28) 



in this region (point B in Fig. 5(b)) and compare it 
to the results in the antiphase and paramagnetic phase 
(point A and C in Fig. 5(b)). The results are plotted in 
Figs. 5(c), 5(d) and 5(e). As can be seen in the plots, 
the spin-spin correlation function at point A (Fig. 5(c)) 
exhibits perfect long-range order which is characteristic 
of the antiphase. At point C in the paramagnetic phase, 
the spin-spin correlation function (Fig. 5(e)) decays ex- 
ponentially with spin separation. At point B in the green 
area in the phase diagram, the spin-spin correlation func- 
tion (Fig. 5(d)) appears to decay algebraically which is 
indicative of the floating phase. 



V. SIMULATING THE ANNNI MODEL IN THE 
SIX-JUNCTION ARRAY SYSTEM 

The achievable scale of an experimental simulation sys- 
tem is limited, by both decoherence and the requirement 
for the two-level approximation to be valid (see appendix 
A). In order to study the feasibility of our Josephson cir- 
cuit system for simulation of frustrated physics, we ex- 
amine a small system to see how close its ground state is 
to certain phases in Fig. 5(b). This information will help 
us determine if it is possible to study the essential char- 
acteristics of a frustrated spin system using a quantum 
simulator of limited size. 

Taking a 6-site Josephson junction array as an exam- 
ple, we obtain the ground state \ip)g of the corresponding 

= 6 ANNNI model using exact diagonalization. Such 
a short chain is insufficient to exhibit the characteris- 
tics of the floating phase, therefore we will focus on the 
ferromagnetic, paramagnetic and antiphase phases. We 
calculate the probabilities oi\ip)g being the ferromagnetic 
and paramagnetic states. 
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P{PM) 
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The results are plotted in Fig. [5] for different values of B 
and A. 

We can see in Fig. [5] that there exists a clear junction 
point A = 0.5, which agrees well with the result in Fig. 
5(b). If we associate the high probability regime with the 
corresponding phase, we can see that the phase regimes 
are nearly the same as well. This indicates that the 6-site 
example can already reveal some essential properties of 
the ANNNI model. 
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FIG. 6: (Color online) (a) The probability of the ground state 
\tp}g of the 6-site Josephson junction array being the ferro- 
magnetic state, (b) The probability of \^)g being the param- 
agnetic state. 
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FIG. 7: The fidelity of the system state to the state of the 
corresponding ANNNI model, for different number of dynam- 
ical decoupling control sequences. 



The results in Fig. |6] are based on exact diagonaliza- 
tion of the AT = 6 ANNNI Hamihonian. In order to 
obtain the ANNNI model from the Hamiltonian of the 
Josephson junction array, dynamical decoupling pulse se- 
quences need to be applied to the qubits as shown in Fig. 
[3l To study the error of the pulse engineered ANNNI 
model, we evaluate the evolution of the Josephson junc- 
tion array system in a time of T = tt under the con- 
trol pulses, and compare it to that of a strict ANNNI 
model in the same amount of time. The initial state 
of the system is set to the maximal superposition state 
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FIG. 8: (Color online) (a) Phase Transition process of the 
ANNNI model with A = 0.4. (b) Phase transition process of 
the Josephson junction array system with and without the 
dynamical-decoupling control sequence. The black dashed 
curve is the result without dynamical-decoupling control se- 
quences. The solid blue and dotted red curves are the results 
with 1 and 4 control sequences in unit time. 



l*o> = 0(1 t>» + I i>0/V2. We divide the total time 

i 

T into m identical time interval: T — mdt. In each in- 
terval 6t, the dynamical decoupling sequences shown in 
Fig. 13] are applied. In Fig. [71 the fidelities of the Joseph- 
son junction array system state compared to the state 
of the corresponding ANNNI model is plotted. We see 
that the fidelities increase with the number m of decou- 
pling sequences, and they already reach 95% within a few 
sequences. 

We can further evaluate the effectiveness of our pulse 
control scheme by comparing how the ground state 
changes with the external field B ina. strict 6-site ANNNI 
model and in a pulse sequence controlled 6-site Joseph- 
son junction array. We simulate this process by adiabat- 
ically changing B in time and calculating the probability 
of the system ground state being in the ferromagnetic 
phase. Here, we set A = 0.4. In Fig. 8(a), we calculate 
the state evolution of a 6-site ANNNI system initially 
in the ferromagnetic state | H ■ ■ ■)z when B — 0. We 
gradually increase the magnetic field B from with veloc- 
ity V = 0.002, and plot the probability that the system 
remains in the original ferromagnetic state at different 
values of B. We can see that when B becomes large, the 
system has deviated from the ferromagnetic state, indi- 
cating that it has changed to a different phase. In Fig. 
8(b), we carry out the same study on the 6-site Josepshon 
junction array under the control pulse scheme in Fig. [31 
As can be seen, the probability of the system remains 
in the ferromagnetic state is almost identical to that in 
the corresponding ANNNI model when the number of 
control sequences is sufficient. These studies show that 
our control pulse based scheme can be used to accurately 
simulate the frustrated ANNNI model. 



VI. CONCLUSION 

In conclusion, we have shown how to simulate frus- 
trated spin models using strongly coupled Josephson 
junction array. We find that the system Hamiltonian 
is exactly solvable beyond the small coupling limit, and 
we design a dynamical decoupling scheme to engineer 
the Hamiltonian for quantum simulation of the ANNNI 
model. We calculate the phase diagram of the system 
numerically using the TEBD method, and demonstrate 
that our control pulse based scheme can be used to sim- 
ulate the corresponding ANNNI model accurately. 
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Appendix A: Conditions for two-level approximation 

In deriving the spin system Hamiltonian in Eq. (|13p . 
we used the two-level approximation for the Josephson 
qubits which kept only the rii = 0, 1 states. We study 
the conditions for the two-level approximation to remain 
valid in this Appendix. 

The Hamiltonian for the Josephson junction array sys- 
tem considering contributions from all charge states reads 

1 ^ 

Hn - -^(n,- 1/2)2 

2—1 

+EEw'("^- V2)(n^+,-i/2) 

i j=l 

-BY,Y.^\n),{n + l\ + h.c.), (Al) 

i n 

where terms in the first line are the on-site charging en- 
ergies of the Josephson qubits, terms in the second line 
are coupling energies between qubits, and terms in the 
third line are the Josephson tunneling energies. When 
the effective magnetic field B is nonzero, the Josephson 
tunneling energies can potentially cause leakage out of 
the rii = 0, 1 states and invalidate the two-level approxi- 
mation. 

We take the ferromagnetic phase of the ANNNI model 
as an example to estimate the probability for the c^ubits 
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in the system to escape the = 0, 1 states. Initially, as- 
sume B = and the system is in the ferromagnetic state 
|5')g = |0101 . . .) (recall that there has been a canoni- 
cal transformation applied on the even sites.). When B 
increases , the qubits can make transitions out of the 
Ui = 0,1 states. We examine the system states that re- 
sult when one of the qubits originally in the n = 1 state 
changes to the n = 2 state, because they are closest to 
|\E')g in energy among all states that violate the two-level 
approximation. We denote such states |4')n- Their en- 
ergies differ from that of |^')g by AE = 2{1/X-1 + A) 
according to Eq. (|A1[) . 

According to the first order perturbation theory, the 
ground state with a nonzero magnetic filed is 



I*) 



Eg — En 



B 



= i^>s + ^(s:j^)„). 



(A2) 



With the form of \'^')g, the total escaping probability 
out of the rii = 0, 1 Hilbert subspace can be estimated to 
be 



B .2iV 



(A3) 



where the factor of N is due to the translational symme- 
try. 



From the Eq. (|A3|1 . we find that the escaping probabil- 
ity is proportional to N. Therefore the allowable system 
size TV is limited by the tolerable escaping probability and 
the amplitude of the magnetic field. For example, if an es- 
caping probability of 5% is acceptable, and B = 0.2, the 
allowable size of the Josephson junction array is about 
N = 10. 
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